﻿using System;
using System.Collections.Generic;
using System.Linq;

namespace MathService.Temp
{
    /// <summary>
    /// Данный класс содержит реализацию симплекс-метода.
    /// </summary>
    /// <remarks>
    /// Изначально был написан на C++, источник: http://marioslapseofreason.blogspot.com/2010/12/linear-programming-simplex-algorithm.html
    /// Портирование под .NET и добавление проверки начальных условий: Владимир Панченко
    /// </remarks>
    public static class Simplex
    {
        /// <summary>
        /// Проверяет корректность исходных данных.
        /// </summary>
        /// <param name="b">Вектор со значениями правой части неравенств.</param>
        /// <param name="MaxFunction"></param>
        private static void CheckSimplex(double[] MaxFunction, double[,] A, double[] b)
        {
            if(MaxFunction.Length != A.GetLength(1))
                throw new ArgumentException("Количество переменных в целевой функции и в матрице коэффициентов не совпадает.");

            if(A.GetLength(0) != b.Length)
                throw new ArgumentException("Количество ограничений и количество значений правой части неравенств не совпадает.");
            for(var i = 0; i < b.Length; i++)
                if(b[i] <= 0)
                    throw new ArgumentOutOfRangeException(nameof(b), "Не выполнено условие: b > 0.");
        }

        /// <param name="maxFunction">Vector with the values of the coefficients for each variable.
        /// Example: { 2 , 1.5 , 1 } for 2X + (3/2)Y + Z</param>
        /// <param name="A">Matrix with the coefficient of the left sides of all the inequalities.</param>
        /// <param name="b">Vector with the values of the right side of the inequalities.</param>
        /// <returns>Linear program in canonical form</returns>
        private static double[][] SetSimplex(double[] maxFunction, double[,] A, double[] b)
        {
            //double[,] simplex;
            var simplex = new List<double[]>();

            var numVariables = maxFunction.Length;
            var numEquations = A.GetLength(0);
            var numCols = numVariables + numEquations + 1 + 1;


            for(var j = 0; j < numEquations; j++)
            {
                var row = new double[numCols];
                for(var i = 0; i < numVariables; i++)
                    row[i] = A[j, i];
                row[numVariables + j] = 1;
                row[numCols - 1] = b[j];

                simplex.Add(row);
            }


            var lastRow = new double[numCols];
            for(var i = 0; i < numVariables; i++)
                lastRow[i] = 0 - maxFunction[i];
            lastRow[numVariables + numEquations] = 1;
            simplex.Add(lastRow);

            return simplex.ToArray();
        }

        private static bool GetPivots(double[][] simplex, ref int pivotCol, ref int pivotRow, ref bool noSolution)
        {
            var numRows = simplex.GetLength(0);
            var numCols = simplex[0].Length;
            noSolution = false;

            var min = 0.0;
            for(var i = 0; i < numCols - 2; i++)
            {
                var value = simplex[numRows - 1][i];
                if(value >= min) continue;
                pivotCol = i;
                min = value;
            }


            if(min.Equals(0))
                return false;


            var minRatio = 0.0;
            var first = true;
            for(var i = 0; i < numRows - 1; i++)
            {
                var value = simplex[i][pivotCol];
                if(value <= 0.0) continue;
                var colValue = simplex[i][numCols - 1];
                var ratio = colValue / value;
                if(!first && ratio >= minRatio || ratio < 0.0) continue;
                minRatio = ratio;
                pivotRow = i;
                first = false;
            }

            noSolution = first;
            return !noSolution;
        }

        /// <summary>
        /// Максимизирует целевую функцию с учётом ограничений с помощью симплекс-метода.
        /// </summary>
        /// <param name="maxFunction">Вектор коэффицентов целевой функции.</param>
        /// <param name="A">Матрица коэффициентов при параметрах в левой части неравенств.</param>
        /// <param name="b">Вектор со значениями правой части неравенств.</param>
        /// <param name="max">Максимальное значение функции с учётом ограничений.</param>
        /// <returns>Вектор значений параметров.</returns>
        public static double[] DoSimplex(double[] maxFunction, double[,] A, double[] b, ref double max)
        {
            CheckSimplex(maxFunction, A, b);

            var simplex = SetSimplex(maxFunction, A, b);

            int pivotCol = 0, pivotRow = 0;
            var numRows = simplex.Length;
            var numCols = simplex[0].Length;


            var noSolution = false;
            while(GetPivots(simplex, ref pivotCol, ref pivotRow, ref noSolution))
            {
                var pivot = simplex[pivotRow][pivotCol];
                for(var i = 0; i < numCols; i++)
                    simplex[pivotRow][i] /= pivot;
                for(var j = 0; j < numRows; j++)
                {
                    if(j == pivotRow) continue;
                    var ratio = -1 * simplex[j][pivotCol];
                    for(var i = 0; i < numCols; i++)
                        simplex[j][i] = simplex[pivotRow][i] * ratio + simplex[j][i];
                }
            }


            if(noSolution)
            {
                var vec = new double[0];
                return vec;
            }

            //optimo!!!
            max = simplex[numRows - 1][numCols - 1];
            var numVariables = numCols - numRows - 1;
            var x = Enumerable.Repeat<double>(0, numVariables).ToArray();

            for(var iCol = 0; iCol < numVariables; iCol++)
            {
                var isUnit = true;
                var first = true;
                var value = 0.0;
                for(var j = 0; j < numRows; j++)
                {
                    if(Equals(simplex[j][iCol], 1.0) && first)
                    {
                        first = false;
                        value = simplex[j][numCols - 1];
                    }
                    else if(!Equals(simplex[j][iCol], 0.0))
                        isUnit = false;
                }
                x[iCol] = isUnit && !first ? value : 0.0;
            }


            return x;
        }
    }
}